rm(list=ls())
library(plyr)
library(foreign)
library(reshape2)
library(readstata13)
#function for standard error of a mean
mean.se <- function (x, na.rm) {stderr <- sd(x, na.rm = TRUE) 
	return(stderr/sqrt(sum(!is.na(x))))}
as.numeric.factor <- function(x) {as.numeric(levels(x))[x]}
options(digits = 3)
options(scipen = 100)

#read in poll data 
megapoll <- read.dta("gay_marriage_megapoll_updated.dta")

#READ IN STATE ESTIMATES
mrp.data <- read.csv('state_preds_means_quantiles_from_STAN.csv')
mrp.data$X <- NULL
mrp.data$state.op <- mrp.data$mrp.mean
with(mrp.data,mean(state.op[state=="NY" & year==1993]))
n.years <-  length(unique( mrp.data$year))
un.year <- unique( mrp.data$year)
un.state <- unique( mrp.data$state)
#plot over time --rder by lowest in 2015
states.sorted  <- arrange(subset(mrp.data, year==2015), mrp.mean)$state

#this is similar to Figure 5, but only contains public opinion, and not policy/court decisions
#pdf("Figure_5.pdf", height=10, width=9)

axis.text <- .9
#par(mfrow=c(10,5),mar=c(1,2,1,1))
layout( matrix(c(1:55), nrow=11, ncol=5, byrow=TRUE) , 
	heights=c(rep(1,10),.3))
par(mar=c(.5,1.5,.5,.7))
for(i in 1:50){
	ok <- mrp.data$state==states.sorted[i]
	plot(mrp.data$year[ok], mrp.data$state.op[ok], type="n", ylim=c(0,80),axes=FALSE, xlab="", ylab="")
	polygon(c(rev(mrp.data$year[ok]), mrp.data$year[ok]),
			 c(rev(mrp.data$mrp.trend.lower[ok]), mrp.data$mrp.trend.upper[ok]),col = 'grey80', border = NA)
	points(mrp.data$year[ok], mrp.data$state.op[ok], type="l", col="blue")
  mtext(states.sorted[i],1,line=-1, adj=.9)
#	axis(1, mgp=c(2,.5,0), at=c(1996,seq(2000,2015,5)))
	axis(2, mgp=c(2,.5,0),las=1,at=seq(0,100,25))
	box()
	abline(h=50, lty=2)
	}
par(mar=c(0,1.5,0,.7))
for (i in 1:5){
	plot(mrp.data$year[ok], mrp.data$state.op[ok], type="n", xlim=range(mrp.data$year), ylim=c(0,80),axes=FALSE, xlab="")
	axis(1, mgp=c(2,.3,0), at =seq(1993,2015,1), line=-2, labels=FALSE)
	axis(1, mgp=c(2,1,0), at =seq(1995,2015,5),las=1,line=-2, tcl=-1.2, cex.axis=axis.text)
	}
	
#dev.off()

#now compare disagg 
state.yearly.disagg <- ddply(megapoll,.(year,state), summarise, 
	mean.state.op.disagg = 100*mean(favor_marriage_comb),
	mean.se = 100*mean.se(favor_marriage_comb),
	mean.lower =	mean.state.op.disagg -1.96*mean.se,
	mean.upper =	mean.state.op.disagg +1.96*mean.se,
	n.obs = length(favor_marriage_comb))
state.yearly.disagg$mean.lower <- ifelse(state.yearly.disagg$mean.lower <0,0,state.yearly.disagg$mean.lower )
state.yearly.disagg$mean.upper <- ifelse(state.yearly.disagg$mean.lower >100,100,state.yearly.disagg$mean.upper)


#PLOT CORRELATION BETWEEN MRP AND DISAGG
#compare correlation
#how many state-combose have fewer than 20 respondents?
temp <- join(mrp.data , state.yearly.disagg)
mean(temp$n.obs < 20, na.rm=T)
mean(temp$n.obs < 10, na.rm=T)

pdf("Figure_A1.pdf", height=5, width=10)

par(mfrow=c(1,2), mar=c(3,3,.5,.5), pty="s")
text.size <-1.2
axis.text <- 1.2
rho.size <- 1.4

plot(temp$state.op, temp$mean.state.op.disagg, xlim=c(0,100), ylim=c(0,100), xaxs="r", yaxs="r", axes=FALSE, pch=21, xlab="", ylab="")
axis(1, mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)
abline(0,1, lty=1)
mtext("MRP estimates of support for gay marriage", 1, cex=text.size, line=1.5)
mtext("Disaggreation estimates", 2, cex=text.size,  line=2, srt=90)
rho <- cor(temp$state.op, temp$mean.state.op.disagg, use="complete.obs")
text(10, 82, expression(paste(bold(rho), "= ")), cex=rho.size)
text(12, 83, round(rho,2), cex=rho.size, adj=0)
text(80,5, "A) All states", cex=text.size, font=4)
box()

#now look at observations with at least 50 resondents) #states with the largest number of respondents
#this is to pull out 10 biggest states
# foo <- tapply(megapoll$restrict, megapoll$state, length)
# foo <-names(sort(foo))[41:50]
# ok <- temp$state %in% foo
ok <- temp$n.obs>50
sort(table(temp$state[ok]))

plot(temp$state.op[ok], temp$mean.state.op.disagg[ok], xlim=c(0,100), ylim=c(0,100), xaxs="r", yaxs="r", axes=FALSE, pch=21, xlab="", ylab="")
axis(1, mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)
abline(0,1, lty=1)
mtext("MRP estimates of support for gay marriage", 1, cex=text.size, line=1.5)
mtext("Disaggreation estimates", 2, cex=text.size,  line=2, srt=90)
rho <- cor(temp$state.op[ok], temp$mean.state.op.disagg[ok], use="complete.obs")
text(20, 82, expression(paste(bold(rho), "= ")), cex=rho.size)
text(22, 83, round(rho,2), cex=rho.size, adj=0)
text(80,5, "B) States with >\n50 respondents", cex=text.size, font=4)
box()

dev.off()

#validate against gay marriage referendum
mrp.data$stateyear <- as.factor(paste(mrp.data$state, mrp.data$year))
state.yearly.disagg$stateyear <- as.factor(paste(state.yearly.disagg$state, state.yearly.disagg$year))

ref.data <- read.csv('state_referendum_results.csv')
#how many ref/states?
dim(ref.data)
length(unique(ref.data$state))
range(ref.data$year)
ref.data$stateyear <- as.factor(paste(ref.data$state, ref.data$year))
ref.data <- merge(mrp.data, ref.data)
ref.data$diff <- abs(ref.data$oppose_ban-ref.data$state.op)
ref.data <- arrange(ref.data, -diff)
#biggest discrepancie occur among states in which civil unions bans were also included
#add disagg
ref.data <- merge (ref.data, state.yearly.disagg)

pdf("Figure_A2.pdf", height=7, width=12)

max.margin<-70
rho.size <- 1.6
states.size <- 1
layout(rbind(c(1:4)), widths=c(.1,1,1,1))
par(mar=c(0,0,0,0))
plot(ref.data$mean.state.op.disagg,ref.data$oppose_ban,xlim=c(10,max.margin), ylim=c(10,max.margin), type="n", xaxs="r", yaxs="r", axes=FALSE, pch=19, xlab="", ylab="")
mtext("Opposition to referenda banning gay marriage", 2, cex=text.size,line=-2,  srt=90,xpd=NA)

par(mar=c(3,2,.5,1),  pty="s")
plot(ref.data$mean.state.op.disagg,ref.data$oppose_ban,xlim=c(10,max.margin), ylim=c(10,max.margin), type="p", xaxs="r", yaxs="r", axes=FALSE, pch=19, xlab="", ylab="")
text(ref.data$mean.state.op.disagg+1, ref.data$oppose_ban,  ref.data$stateyear, cex=states.size, adj=0)
axis(1, mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)
abline(0,1, lty=1)
rho <- cor(ref.data$oppose_ban, ref.data$mean.state.op.disagg, use="complete.obs")
text(20, 62, expression(paste(bold(rho), "= ")), cex=rho.size)
text(21, 62.5, round(rho,2), cex=rho.size, adj=0)
box()
mtext("Disaggregated estimates\nof support for gay marriage", 1, cex=text.size, line=4)
#mtext("Opposition to referenda banning gay marriage", 2, cex=text.size, line=2, #srt=90)
mtext("A) Disaggregation", 3, cex=2, font=2, line=1)

max.margin <- 70
plot(ref.data$state.op,ref.data$oppose_ban,xlim=c(10,max.margin), ylim=c(10,max.margin), type="p", xaxs="r", yaxs="r", axes=FALSE, pch=19, xlab="", ylab="")
text(ref.data$state.op+1, ref.data$oppose_ban,  ref.data$stateyear, cex=states.size, adj=0)
axis(1, mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)
abline(0,1, lty=1)
rho <- cor(ref.data$oppose_ban, ref.data$state.op, use="complete.obs")
text(20, 62, expression(paste(bold(rho), "= ")), cex=rho.size)
text(21, 62.5, round(rho,2), cex=rho.size, adj=0)
box()
mtext("MRP estimates of support\nfor gay marriage", 1, cex=text.size, line=4)
#mtext("Opposition to referenda banning gay marriage", 2, cex=text.size,  line=2, srt=90)
mtext("B) MRP (all referenda)", 3, cex=2, font=2, line=1)

#only in same-sex marriage states
ref.data.marr <- subset(ref.data, marr_only==1)
ref.data.marr <- arrange(ref.data.marr, diff)
plot(ref.data.marr$state.op,ref.data.marr$oppose_ban,xlim=c(10,max.margin), ylim=c(10,max.margin), type="p", xaxs="r", yaxs="r", axes=FALSE, pch=19, xlab="", ylab="")
text(ref.data.marr$state.op+1, ref.data.marr$oppose_ban,  ref.data.marr$stateyear,cex=states.size, adj=0)
axis(1, mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, mgp=c(2,.5,0), cex.axis=axis.text, las=1)
abline(0,1, lty=1)
rho <- cor(ref.data.marr$oppose_ban, ref.data.marr$state.op, use="complete.obs")
text(20, 62, expression(paste(bold(rho), "= ")), cex=rho.size)
text(21, 62.5, round(rho,2), cex=rho.size, adj=0)
box()
mtext("MRP estimates of support\nfor gay marriage", 1, cex=text.size, line=4)
#mtext("Opposition to referenda banning gay marriage", 2, cex=text.size,  line=2, srt=90)
mtext("C) MRP (marriage only)", 3, cex=2, font=2, line=1)

dev.off()

